# Load Packages
library(binsreg)

source("code/config.R")
source(paste0(CODE,"exacthat.R"))
source(paste0(CODE,"parameters.R"))
source(paste0(CODE,"binscatter2.R"))

df<-read_csv(paste(DATA, "clean_data4.csv", sep=""))
df$one<-1

# plot change in transit time: Figure A6
df$t_hat<- exp(df$t_post-df$t)
ggplot(df, aes(x=t_post - t))+ geom_density() +ylab("Density") + xlab("Change in Travel Time from DTL (Mins)")

# Function to extract results for table

table <- function(out){
  tab<-matrix(nrow=5, ncol=2)
  tab[1,1:2]<- c(out[[2]],out[[1]])-1 # Welfare
  tab[2,1:2]<- c(out[[15]],out[[14]])-1 # Access to Consumption
  tab[3,1:2]<- c(out[[17]],out[[16]])-1 # Access to Employment
  tab[4,1]<- c(out[[1]])- c(out[[2]]) # Gap in Welfare Impact
  tab[5,1]<- seg(out)-1 # Segregation
 return(tab)
}

# Function to extract gaps
gaps <- function(out){
  return(c(out[[1]]- out[[2]], out[[14]]- out[[15]], out[[16]]- out[[17]]))
}

# Segregation calculation
seg<-function(out_df){
  df$hat_share_commute_N<-out_df[[3]]
  df$hat_share_commute_T<-out_df[[4]]
  df$hat_share_commute_w_N<-out_df[[5]]
  df$hat_share_commute_w_T<-out3[[6]]
  df$hat_share_leisure<-out_df[[7]]
  df$hat_share_leisure_w<-out_df[[8]]
  df$hat_share_resident<-rep(out_df[[9]],each=40)
  df$hat_share_resident_w<-rep(out_df[[10]],each=40)
  df_o<- df %>% select(hat_share_resident,hat_share_resident_w,share_resident,share_resident_w,sub_area_o) %>% distinct()
  sum(abs(df_o$share_resident*df_o$hat_share_resident - df_o$share_resident_w*df_o$hat_share_resident_w))/sum(abs(df_o$share_resident - df_o$share_resident_w))
}

# Estimate counter factual: Table 2

## Main estimates

table(counter_agg(df$t_hat ,df$t_hat)) # Full model
table(counter_c_agg(df$one ,df_pre$t_hat)) # Only commuting

## No spillovers
table(counter(df_pre$t_hat ,df_pre$t_hat)) # Full model
table(counter_c(df_pre$one ,df_pre$t_hat)) # Only commuting

# Decomposition: Table 3

## baseline
r0<-gaps(counter(df_pre$t_hat ,df_pre$t_hat))

## equalize consumption shares
source(paste0(CODE,"parameters.R"))
alpha_N_m<-alpha_N_p
alpha_H_m<-alpha_H_p
alpha_T_m<-alpha_T_p
r1<-gaps(counter(df_pre$t_hat ,df_pre$t_hat))

## equalize travel costs
source(paste0(CODE,"parameters.R"))
kappa_w<-kappa
r2<-gaps(counter(df_pre$t_hat ,df_pre$t_hat))

## equalize travel elasticities
source(paste0(CODE,"parameters.R"))
epilson_l_w<-epilson_l
epilson_s_w<-epilson_s
r3<-gaps(counter(df_pre$t_hat ,df_pre$t_hat))

## equalize travel elasticities
source(paste0(CODE,"parameters.R"))
epilson_R_m<-epilson_R_p
r4<-gaps(counter(df_pre$t_hat ,df_pre$t_hat))

## Generate Table
rbind(r0,r1,r2,r3,r4)


# Extract variables
df_pre$hat_share_commute_N<-out3[[3]]
df_pre$hat_share_commute_T<-out3[[4]]
df_pre$hat_share_commute_w_N<-out3[[5]]
df_pre$hat_share_commute_w_T<-out3[[6]]
df_pre$hat_share_leisure<-out3[[7]]
df_pre$hat_share_leisure_w<-out3[[8]]
df_pre$hat_share_resident<-rep(out3[[9]],each=length(unique(df_pre$sub_area_d)))
df_pre$hat_share_resident_w<-rep(out3[[10]],each=length(unique(df_pre$sub_area_d)))
dest_data<-df_pre %>% group_by(sub_area_d) %>% summarise(shoppers= sum(share_resident*hat_share_resident*share_leisure*hat_share_leisure)/sum(share_resident*share_leisure), shoppers_w= sum(share_resident_w*hat_share_resident_w*share_leisure_w*hat_share_leisure_w)/sum(share_resident_w*share_leisure_w),workersn= sum(share_resident*hat_share_resident*share_commute*hat_share_commute_N*share_N_p + share_resident*hat_share_resident*share_commute*hat_share_commute_T*share_T_p)/sum(share_resident*share_commute), workersn_w= sum(share_resident_w*hat_share_resident_w*share_commute_w*hat_share_commute_w_N*share_N_m + share_resident_w*hat_share_resident_w*share_commute_w*hat_share_commute_w_T*share_T_m)/sum(share_resident_w*share_commute_w))

# Create map of change in Consumption Trips: Figure 1 Panel b
mpsz <- st_read(dsn = paste0(DATA,"master-plan-2014-subzone-boundary-web-shp"), 
                layer = "MP14_SUBZONE_WEB_PL")
DTLStationCoords <- read.csv(paste0(DATA,"DTLStationCoords.csv"), header = TRUE, stringsAsFactors = FALSE)
stations1 <- st_as_sf(DTLStationCoords, coords = c('LON', 'LAT'), crs=4326)
dest_data$PLN_AREA_N<-dest_data$sub_area_d
mpsz<-left_join(mpsz,dest_data)
mpsz$`Percent Change`<-log(mpsz$shoppers)
qtm(mpsz, fill="Percent Change",fill.n = 10) + tm_legend(show=T) +tm_shape(stations1)+tm_bubbles(col = "red", scale = .3) 


## Test predictions ##

df<-read_csv(paste(DATA, "clean_data4.csv", sep=""))
df_post<-read_csv(paste(DATA, "clean_data_post.csv", sep=""))

# Estimate 
df$t_hat<- exp(df$t_post-df$t)
counter_agg(df$t_hat ,df$t_hat)

# Figure A19a
pred_change_leisure<-out[[7]]
change_leisure<-df_post$share_leisure/df$share_leisure
change_leisure[change_leisure==Inf]<-NA
binsreg(change_leisure,pred_change_leisure)
temp<-data_frame(change_leisure=change_leisure,pred_change_leisure=pred_change_leisure)
g<-binscatter2(formula="change_leisure ~ pred_change_leisure", key_var = "pred_change_leisure", data=temp, bins=20, partial=FALSE)
g + xlab("Model: Predicted Change in Consumption Travel") + ylab("Data: Change in Consumption Travel")

# Figure A19b
pred_change_commute<-out[[2]]*df$share_N_p+out[[3]]*df$share_T_p
change_commute<-df_post$share_commute/df$share_commute
binsreg(change_commute,pred_change_commute)
temp<-data_frame(change_commute=change_commute,pred_change_commute=pred_change_commute)
g<-binscatter2(formula="change_commute ~ pred_change_commute", key_var = "pred_change_commute", data=temp, bins=22, partial=FALSE)
g + xlab("Model: Predicted Change in Commuting") + ylab("Data: Change in Commuting")

# Figure A19c
pred_change_res<-out3[[9]]
change_fes<-df_post$share_resident[(1:length(unique(df$sub_area_d)))*length(unique(df$sub_area_d))]/df$share_resident[(1:length(unique(df$sub_area_d)))*length(unique(df$sub_area_d))]
plot(change_fes,pred_change_res)
temp<-data_frame(change_fes=change_fes,pred_change_res=pred_change_res)
g<-binscatter2(formula="change_fes ~ pred_change_res", key_var = "pred_change_res", data=temp, bins=28, partial=FALSE)
g + xlab("Model: Predicted Change in Resident Share") + ylab("Data: Change in Resident Share")
summary(lm(change_fes~pred_change_res))

# Read in food establishments data
neaEatingEstab2015 <- read_dta("data/neaEatingEstab2015.dta")
neaEatingEstab2018 <- read_dta("data/neaEatingEstab2018.dta")
neaEatingEstab2019 <- read_dta("data/neaEatingEstab2019.dta")
neaEatingEstabMerged <- read_dta("data/neaEatingEstabMerged.dta")
neaEatingEstabMerged <- neaEatingEstabMerged %>% select(LICENCE_NUMBER,POSTCODE=POSTAL_CODE,SUBZONE_N,PLN_AREA_N)
neaEatingEstab2015<- neaEatingEstab2015 %>% select(LICENCE_NUMBER,POSTCODE) %>% mutate(y2015=1)
neaEatingEstab2018<- neaEatingEstab2018 %>% select(LICENCE_NUMBER=LIC_NO,POSTCODE) %>% mutate(y2018=1)
neaEatingEstab2019<- neaEatingEstab2019 %>% select(LICENCE_NUMBER,POSTCODE=POSTAL_CODE) %>% mutate(y2019=1)
neaEatingEstab2018$POSTCODE<-as.numeric(neaEatingEstab2018$POSTCODE)

# clean data
alleat<-full_join(neaEatingEstab2015,neaEatingEstab2018)
alleat<-full_join(alleat,neaEatingEstab2019)

alleat <- left_join(alleat,neaEatingEstabMerged)
alleat$y2015[is.na(alleat$y2015)]<-0
alleat$y2018[is.na(alleat$y2018)]<-0
alleat$y2019[is.na(alleat$y2019)]<-0

alleat$enter2018<-(alleat$y2018==1 & alleat$y2015==0)
alleat$enter2019<-(alleat$y2019==1 & alleat$y2018==0)
alleat$exit2018<-(alleat$y2015==1 & alleat$y2018==0)
alleat$exit2019<-(alleat$y2018==1 & alleat$y2019==0)

alleat <-alleat %>% group_by(SUBZONE_N) %>% summarise( y2015= sum(y2015),y2018= sum(y2018), y2019= sum(y2019), gross_entry18=sum(enter2018), gross_entry19=sum(enter2019),gross_exit18=sum(exit2018), gross_exit19=sum(exit2019 )  )
alleat$growth<-alleat$y2018/alleat$y2015
alleat<-alleat %>% filter(SUBZONE_N %in% unique(df_pre$sub_area_o))
pred_change_commercial<-out3[[11]]

# Figure A19d
temp<-data_frame(change_commercial=alleat$growth,pred_change_commercial=pred_change_commercial,entry=alleat$entry_rate, size=df_pre$commercial_quant[(1:40)],estabs=alleat$y2015,name=alleat$PLN_AREA_N)
g<-binscatter2(formula="entry ~ pred_change_commercial", key_var = "pred_change_commercial", data=temp[temp$pred_change_commercial<1.4,], bins=20, partial=FALSE)
g+ xlab("Model: Predicted Change in Non-Tradables Production") + ylab("Food Establishment Entry (2015 - 2018)")



